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Abstract. We consider the primal problem of finding the zeros of the sum of a maximally 
monotone operator with the composition of another maximally monotone operator with 
a linear continuous operator and a corresponding dual problem formulated by means of 
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o 



the inverse operators. A primal-dual splitting algorithm which simultaneously solves the 
two problems in finite-dimensional spaces is presented. The scheme uses at each iteration 
separately the resolvents of the maximally monotone operators involved and it gives rise 
i-^h to a splitting algorithm for finding the zeros of the sum of compositions of maximally 

monotone operators with linear continuous operators. The iterative schemes are used for 
solving nondifferentiable convex optimization problems arising in image processing and in 
location theory. 
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1 Introduction and preliminaries 

y— i In this paper we propose an iterative scheme for solving the inclusion problem 

.£h find x £ X such that G Ax + K*BKx, 

X 

where X and Y are Hilbert spaces, A : X X and B : Y Y are maximally monotone 
operators and K : X — > Y is a linear continuous operator, which makes separately use 
of the resolvents of A and B. The necessity of having such an algorithm is given by the 
fact that the classical splitting algorithms have considerable limitations when employed 
on the inclusion problem under investigation in its whole generality. Indeed, the forward- 
backward algorithm (see [5]) is a valuable option in this sense when B is single- valued and 
cocoercive, while the use of Tseng's algorithm (see [22]) asks for B being single- valued and 
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Lipschitz continuous on a superset of the image of the domain of A through K. On the 
other hand, the Douglas- Rachford algorithm (see [5pl]) asks for the maximal monotonicity 
of A and K*BK and employs the resolvent of the latter, which can be expressed by means 
of the resolvent of B only in some very exceptional situations (see (5j Proposition 23.23]). 

The aim of this article is to overcome this shortcoming by providing a primal- dual 
splitting algorithm for simultaneously solving this inclusion problem and its dual inclusion 
problem in the sense of Attouch-Thera (see [2j[4j[l9]), in the formulation of which the 
resolvents of A and B appear separately. In the case when A and B are subdifferentials 
of proper, convex and lower semicontinuous functions we rediscover as particular case the 
iterative method from [9j. We also show how the provided primal-dual algorithm gives 
rise to a primal-dual iterative method for finding the zeros of the sum of compositions 
of maximally monotone operators with linear continuous operators. The latter will find 
application when solving nondifferentiable convex optimization problems arising in image 
processing and in location theory having in the objective the sum of (more than two) 
compositions of proper, convex and lower semicontinuous functions with linear continuous 
operators. 

For another primal-dual splitting algorithm for simultaneously solving a primal inclu- 
sion problem and its Attouch-Thera-type dual inclusion problem, recently introduced in 
the literature, we refer the reader to [8|]10|. By using a consecrated product space approach, 
this method basically reformulates the primal-dual pair as the problem of finding the zeros 
of the sum of a maximally monotone operator and a monotone and Lipschitz continuous 
operator, which is then solved by making use of the relaxed version of Tseng's algorithm. 

The structure of the paper is the following. The remaining of this section is dedicated 
to some elements of convex analysis and of the theory of maximally monotone operators. 
In Section [2] we motivate and formulate the primal-dual splitting algorithm for solving 
the problem of finding the zeros of the sum of a maximally monotone operator with the 
composition of another maximally monotone operator with a linear continuous operator 
and its dual problem and investigate its convergence properties. In Section [3] we formu- 
late a primal-dual splitting algorithm for the problem of finding the zeros of the sum of 
compositions of maximally monotone operators with linear continuous operators, while in 
Section [4] we employ the two primal-dual schemes for solving several classes of nondiffer- 
entiable convex optimization problems. Finally, we consider applications of the presented 
algorithms in image deblurring and denoising and when solving the Fermat- Weber loca- 
tion problem. For the latter we compare their performances to the ones of some iterative 
schemes recently introduced in the literature. 

In what follows we recall some elements of convex analysis and of the theory of max- 
imally monotone operators in Hilbert spaces and refer the reader in this respect to the 
books (5}[6}[l4|[2l}[2i). 

Let X be a real Hilbert space with inner product (•,•) and associated norm || • || = 
yf (•, •). For a function f : X — > K, where 1 := MU {±00} is the extended real line, we 
denote by dom/ = {x G X : f(x) < +00} its effective domain and say that / is proper if 
dom/ ^ and f{x) > -00 for all x G X. Let /* : X R, /*(«) = sup x&x {{u, x) - f(x)} 
for all u G X, be the conjugate function of /. The subdifferential of / at x G / _1 (M) is 
the set df(x) := {u G X : f(y) > f(x) + (u,y — x) \/y G X}. We take by convention 
df(x) := 0, if x ^ / _1 (M). For every 7 > and every x G X it holds d{ r yf){x) = ^df(x). 
When Y is another Hilbert space and K : X — > Y a linear continuous operator, then 
K* : Y — > X, defined by (K*y,x) = (y,Kx) for all (x,y) G X x Y, denotes the adjoint 
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operator of K. 

Let C C X be a nonempty set. The indicator function of C, <5<7 : A — )• M, is the 
function which takes the value on C and +00 otherwise. The subdifferential of the 
indicator function is the normal cone of C, that is Nc(x) = {u G A : {u,y — x) < 
Vy G C}, if 1 E C and Nc(x) = for x ^ C. If C is a convex set, we denote by 
sqriC := {1 £ C : Ua>oA(C — x) is a closed linear subspace of X} its strong quasi-relative 
interior. The strong quasi-relative interior of C is a superset of the topological interior of 
C, i.e., intC C sqriC (in general this inclusion may be strict). If X is finite-dimensional, 
than sqri C coincides with the relative interior of C, which is the interior of C with respect 
to the affine hull of this set. For more results relating to generalized interiority-type notions 
we refer the reader to (5]-[7}[2l}[24] . 

For an arbitrary set- valued operator A : X 14 X we denote by Gr A = { (x, u) G X x X : 
u G Ax] its graph, by domA = {x G X : Ax 7^ 0} its domain and by A' 1 : X X its 
inverse operator, defined by (it, x) G Gr A -1 if and only if (x, u) G Gr A. We say that A is 
monotone if (x — y, u — v) > for all (x, u), (y, v) G Gr A. A monotone operator A is said 
to be maximally monotone, if there exists no proper monotone extension of the graph of 
A on X x X. Notice that the subdifferential of a proper, convex and lower semicontinuous 
function is a maximally monotone operator (cf. [20] ) . A single- valued linear operator 
A : X — > X is said to be skew, if (x, Ax) = for all x G X. The skew operators 
are maximally monotone and if they are not identical to zero and the dimension of X is 



greater than or equal to 2, then they fail to be subdifferentials (see 21 ). When / : X — > R 
is a proper, convex and lower semicontinuous it holds (9/) _1 = df*. The resolvent of 
A, J A : X ■=!, X, is defined by J A = (Idx+-4)~\ where ld x ■ X -> A,Idx(x) = x 
for all x G AT, is the identity operator on X. Moreover, if A is maximally monotone, 
then J a '■ X — > X is single- valued and maximally monotone (cf. [5j Proposition 23.7 and 
Corollary 23.10]). For an arbitrary 7 > we have (see [SJ Proposition 23.2]) 

p G Jjax if and only if (p, 7 _1 (x — p)) G Gr A 

and (see [5| Proposition 23.18]) 

J^A + 7J7-1A-1 7 _1 Id x = . (1) 

When / : X — > R is a proper, convex and lower semicontinuous function and 7 > 0, for 
every x G A we denote by prox 7 j(x) the proximal point of parameter 7 of / at x, which 
is the unique optimal solution of the optimization problem 

inf (/(y) + ^-||y-x|| 2 l. (2) 



yex { 27 1 

Notice that J^qj = (Idx +7<9/)~ 1 = pvox^t, thus prox 7 j : A — > X is a single- valued 
operator fulfilling the extended Moreau 's decomposition formula 

prox 7/ +7prox (1/7)/ , 07- 1 ld x = ld x ■ (3) 

Let us also recall that the function / : X — )■ K is said to be strongly convex (with modulus 
7 > 0), if / — ?|| • || 2 is a convex function. 

Finally, we notice that for f = Sc, where C C A is a nonempty closed and convex set, 
it holds 

J 7 at c = JAr c = J 9 8 c = (Idx +AC)" 1 = prox^ = P c , (4) 
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where Pc '■ X — > C denotes the projection operator on C (see (5j Example 23.3 and 
Example 23.4]). 

2 A primal-dual splitting algorithm for finding the zeros of 

A + K*BK 

For X and Y real Hilbert spaces, A : X =4 X and B : Y =4 Y maximally monotone 
operators and A' : X 7 a linear and continuous operator we consider the problem of 
finding the pairs (x,y) £ X x Y fulfilling the system of inclusions 

Kx G B~ 1 y and - K*y G Ax. (5) 

If (x,y) fulfills g), then x is a solution of the primal inclusion problem 

find x e X such that £ Ax + K*BKx (6) 

and y is a solution of its dual inclusion problem in the sense of Attouch-Thera 

find y G Y such that G B~ x y - KA- l (-K*)y. (7) 

On the other hand, if x G X is a solution of the problem (|6]), then there exists a solution y 
of Q such that (x, y) fulfills ([5]), while, if y € Y is a solution of the problem 0, then there 
exists a solution x of ([6]) such that (x, y) fulfills We refer the reader to [TJ[2j|4j|8 13 , 19 



for more algorithmic and theoretical aspects concerning the primal-dual pair of inclusion 
problems ([6])-([7]). 

For all a, r > it holds 

(x, y) is a solution of Q 43- y + oKx G (Idy -I-ct-B -1 )?/ and x — rK*y G (Idx +t^4)x 
^ V = JaB- 1 {y + ctKx) and x G J t a{x — rK*y). 

(8) 

The above equivalences motivate the following algorithm for solving §5§. 
Algorithm 1 

Initialization: Choose a, r > such that crr||ET|| 2 < 1 and (x°,y°) G X x Y. 
Set x° := x°. 



For n > set: y 

x n+l 



n+l ~ J aB -i(y n + aKx n ) 
J rA (x n - TK*y n+1 ) 

2x n+l _ x n 



Theorem 2 Assume that the system of inclusions ^ has a solution (x, y) G X x Y and 
let (x n ,x n ,y n ) n >o be the sequence generated by Algorithm^ The following statements are 
true: 

(i) For any n > it holds 

- -ii 2 + (1 _ aTm f)\\y n -y\\ 2 < ii" - g n 2 + n^-fli 2 , (9 ) 



2r v " " ' 2a ~ 2t 2<t 

i/ms i/ie sequence (x n ,y n ) n >o is bounded. 

(ii) If X and Y are finite-dimensional, then the sequence (x n ,y n ) n >o converges to a 
solution of the system of inclusions 
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Proof, (i) For any n > the iterations of Algorithm [T] yield that 



x n+1 , - rK*y n+1 - x n+1 ) ) G Gri, 



hence the monotonicity of A implies 



1 



< ( x n+1 - x, -{x n - x n+l ) - K*y n+l + K*y 



Similarly, for any n > we have 



y n+ \ -{y n + oKx n - y" +1 ) ) € Gr B~\ 
a 



thus 



< ( Kx n + -{y n - y n+L ) - Kx, y n+l - y > . 



On the other hand, for any n > we have that 



\ x ^-x\\ 2 + (x n+1 



1 



x, -{x n - x n+1 ) - K*y n+l + K*y 

T 



(10) 

(11) 

(12) 
(13) 



x n+1 - x, \x n - x) + ( 1 - - ) {x n+1 - x) - K*y n+1 - K /) 



-(x n+1 -x,x n -x) + ( 1 - 1 ) \\x n+1 -xf + (K(x n+i -x),-y n+1 +y), 



hence 



\x n+i -m 2 + (x r + l 



x, -{x n - x n+1 ) - K*y n+1 + K*y 

T 



- {x n+l -x,x n -x) + (K(x n+1 - x), -y n+l + y) = 

T 

- x n f + - x\\ 2 + - x n f + (K(x n+1 - x), -y n+1 + y), 

2r 2r 2t 



where, for deriving the last formula, we use the identity 

(a, b)- 

Consequently, for any n > it holds 



^ll«-fe|| 2 + ^l|a|| 2 + ^||6|| 2 Va,6GX 



-*-||x n+1 - x\\ 2 + (x n+1 - x, -{x n - x n+1 ) - K*y n+l + K*y 
2r \ t 



~\\x n+1 - x n f + ^-\\x n - xf + {K(x n+1 - x), -y n+l + y). 



(14) 



Thus, by combining (11) and (14), we get for any n > 



-*-||x n+1 - x\\ 2 < ~^-\\x n+1 - x n \\ 2 + - x\\ 2 + (K(x n+1 - x), -y n+l + V). (15) 

2r 2t 2t 



By proceeding in analogous manner we obtain the following estimate for any n > 

1 



,n+l 



n n+l\ 



Kx, y 



n+l 



y 



1 (y n - y) + ( 1 - - ) (y n+1 - y ) + Kx n - Kx, y n+s 



a 



{y n_fr y n+l_~ )+ , 1 



a 



a 



J n+l -y\\ 2 + (K(x n -x),y n+l -y), 



hence 



1 l|y n+1 - VW 2 + (Kx n + -{y n - y n+1 ) - Kx, y n+1 



1 



a 



(y n - y, y n+ ' - y ) + (K(x n - x),y n+1 - y) 



-^\\y n+l -y n t + Y a \\y n+1 

From here we obtain for any n > 



+ 



5' 



y -f + (K(x n -x),y n+1 -y). 



1 

2a 1 



\W n+1 - yf + {Kx n + -(y n - - Kx,y 



n „,n+l\ 



..n+l 



a 



~\\y n+1 -y n \\ 2 + —\\y n 



(16) 



+ {K(x n -x),y n+1 -y). 



and, thus, by combining (13) and (16), it follows 

^h n+1 -yf< -^\\y n+l - y n \? + ^h n 



+ {K{x n -x),y n+l -y). (17) 



Summing up the inequalities (15) and (17) and taking into account the definition of 
xf 1 we obtain for any n > 



^\\x n+1 -^ 2 + ^\\y n+1 -y\\ 2 < 



2t 



x" - x r + 



si 



I 2 - — \\x n+1 - x n \\ 2 - —\\y n+1 - y n f+ 
i 2t \\ ii 2(j \\y y ii -r 



(18) 



(K(x n+1 + x^ 1 - 2x n ), -y n+1 + y) 



where x 1 := x°. Let us evaluate now the last term in relation (18). For any n > it 
holds 

(K(x n+1 + x 11 - 1 - 2x n ), -y n+1 + y) = 
(K(x n+1 - x n ), -y n+1 + y) + (K(x n - x n - 1 ),y n -y) + (K(x n - Z" 1 ), y n+1 - y n ) < 
-(K(x n+1 - x n ),y n+1 ~y) + (K(x n - x n ~ 1 ),y n -y) + 

x n ),y 



x n -s n - 1 ||||y n+1 -i/ n ||< 



-(K(x n+1 - x n ),y n+1 -y) + (K(x n - x^ 1 ), y n - y)- 

V°T\\ K W \\ x n _ x n-l\\2 + V^ll^ll || y n+1 _ „.n||2 



2r 



2a 



y 



(19) 
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From (18) and (19) we obtain for any n > the following estimation 



— ||x n+1 - x\\ 2 + — \\y n+1 - y\\ 2 < —\\x n - xf + — \\y n - yf- 
2r" 11 2a uy " ~ 2r" 11 2a " " 



1 



,1^+1 _ x n ||2 _ 1 yn+1 _ y n {l 2 _ {K ^n+X _ ^ y n+l _ ^ + _ x n-l^ y n _ y) + 

2r 2a 



^T\\ K h \ x n X n-U2 



2a 



thus 



|x n+1 -x|| 2 \\y n+1 



2t 



2a 



I 2 \\x n — xll 2 
- 2t 



-+ 



(-1 + - " n _ " + y/*F\\K\\ 



2a 



(20) 



2a 2t v " " 2r 

<^(x™ +1 - x n ),y n+1 -y) + (K(x n - x n ~ 1 ),y n - y). 



Let be an arbitrary N 6 N, iV > 2. Summing up the inequalities in ( |20| from n = to 
N — 1 we obtain 



,N 



2t 



+ 



2a 



< 



2t 



+ 



2a 



-+ 



N II n n—l\\2 

(-1 + v^n*ii) E + (- 1 + v^ii^id E 

n=l 



' ' X^ X^ "^11^ 



2t 



\ X N _ X N-1\\2 
2r 



n=l 

(K(x N -x N -'),y N -y). 



(21) 



By combining (21) with 



(K(x N -x N - 1 ),y N -y)< 



\x — x 



2t 



+ 



2a 



we get 



Ix^-xH 2 Hy^-yll 2 < ||x°-£|| 2 



2t 



2a 



N || n n— 1 II 2 

E ~2a +(-1+^11^11) E 



2t 2<t 

iV_1 ll„n ^,71.-1 1|2 



+ 



n=l 



X — X 



n=l 



2r 



+ 



<rr\\K\\< 
2a 



\\y N -y\\ 2 



or, equivalently, 



- -ill - . , ,-„ |2n \\y N - y\\ 2 



2t 



+ (I - ar\\K\\y- 



2a 



+ 



N || 7 ,™ _ n ri ~ l \\ 2 

(i - ^m\\k\\) e lly I " + (i - v^\\k\\) e 

n=l n=l 



iV ~ 1 "x n — x n_1 |l 2 



2t 



< (22) 



|x°-x|| 2 + ||y°-y|| 2 



2t 2<t 

By taking into account that ctt|| J^T|| 2 < 1 (22) yields Q, hence (x n ,y n ) n >o is bounded. 



(ii) According to (i), {x n ,y n ) n >o has a subsequence (x Uk , y nk )k>o which converges to 
an element (x*,y*) G X x Y as k — > +00. From (10) and (12) and using that, due to the 
maximal monotonicity of A and B, Gr A and Gr B are closed sets, it follows that (x* , y*) 
is a solution of the system of inclusions ([5]). On the other hand, from (22) we obtain that 
lim n _, +00 (x™ - x"" 1 ) = lim n _, +00 (y™ - y™" 1 ) = 0. 

Further, let be k > and iV E N, jV > n^. Summing up the inequalities in (20), for 
(x, y) := (x* , y*), from n = to N — 1 we obtain 



\x N -x*\\ 2 
27 



„W"_,.*||2 ^ 



2cr 



n fc -ll|2 



jV-1 



2r 



+ (1-^1*11) E 



x — x 



n=nj,+l 
n— 1 112 



\y n -y n ~ l f 

2a 



n=n k 



2t 



+ 



\ X N _ x N-lp 

27 



+ (K(x N - x N ~ 1 ),y N - y*) - (K(x n « - x^ 1 ),^ - y*) 



< 



\ x n k _ x * 112 



2r 



+ 



\y nk -y*\\ 2 

2a 



which yields 



| X N _ x *||2 

27 



1 y V \\ ^ 11 7^1111 jv JV-i 11 11 jv *n 1 
2cr ~ 11 



|x nfc -x*|| 2 l|y nfc -w*|| 2 



+ 



+ 



\x nk — x nfe_1 II 2 



+ (K(x nk -x n "- 1 ),y nk -y* 



2t 2a 2t 

Consequently, by using also the boundedness of (x™,y n ) n >o, for any k > it holds 



lim sup 

jv^+oo 



| X JV_ X *I|2 



2r 



+ 



\y N -y*\\ 2 

2a 



< 



| x ™fc_ x *|| 2 ||y"fc_ y *||2 \\x Uk _ x n k -l\\2 



+ 



2r 2cr 2r 

We finally let k converge to +00, which yields 



+ (K(x nk -x nk - 1 ),y nk -y* 



lim sup 

jv^+oo 



\x N -x*\\ 2 
2t 



\y N -y*\? 

2a 



and, further, limjv->+oo x = x* and lini7v^+oo y = y 



Remark 3 The characterization of the solution of the system of inclusions ([5]) given in 
^ motivates the following iterative scheme 
Initialization: Choose a, r > and (x°,y°) G X x Y. 
For n > set: y n+1 := J aB -i (y n + aKx n ) 



„n+l 



:= J rA (x n - rK*y 



*„.ra+l\ 



as well, which is nothing else than an Arrow-Hurwicz-Uzawa-type algorithm (see [3]) 
designed for the inclusion problem Q. 
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We close this section by discussing another modality of investigating the system of 
inclusions ^ by employing some ideas considered in (8j[l0j. To this end we define the 
operators M: XxY^XxY, M(x,y) = (Ax,B~ l y), and S : X x Y -)• X x Y , 
S(x,y) = (K*y,—Kx). The operator M is maximally monotone, since A and B are 
maximally monotone, while S is maximally monotone, since it is a skew linear operator. 
Then (x,y) E X x Y is a solution of the system of inclusions §5§ if and only if it solves 
the inclusion problem 

find (x, y) G X x Y such that (0, 0) G S(x, y) + M(x, y). (23) 

Applying Algorithm [l] to the problem (23) with starting point (x°, y°, u°, v°) G X x 
Y" x X x Y", (x°,y°) = (x°,y°) and cr, r > gives rise for any n > to the following 
iterations: 

J aM -,[(u n ,v n ) + a(x n ,y n )] 



(u n+1 ,v n+1 ) 
(x n+1 ,y n+1 ) 
(x n+1 ,y n+1 ) 



J rS [(x n ,y n )-r(u 



n+l v n+l 



Since 



2(x n+1 ,y n+1 ) - (x n ,y n ). 

JaM- 1 = JaA' 1 x J<tB 



For n > set: n n+1 

y" +1 



and (see [8j Proposition 2.7]) 

Jrs{x,y) = ((ld x +t 2 K*K)-\x - rK*y), (ld Y +r 2 KK*)- l (y + rKx)) V(x,y) G X x y, 
this yields the following algorithm: 
Algorithm 4 

Initialization: Choose cr, r > such that ar < 1 and (u°,v°) £ X xY. 

Set (x°,y°) := (x°,y°). 

J aA -i(u n + ax") 
JaB(v n + ay n ) 

(Id x +t 2 K*K)- 1 [x n - ru n+1 - rK*(y n - Tv n+1 )] 
(Idy +t 2 KK*)- 1 [y n - rv n+1 + rK(x n - ru n+1 )} 

y • - ■= 2y n+1 - y n 

The following convergence statement is a consequence of Theorem [2} 

Theorem 5 Assume that X and Y are finite- dimensional spaces and that the system 
of inclusions ^ is solvable. Then the sequence (x n ,y n ) n >o generated in Algorithm^ 
converges to (x*,y*), a solution of the system of inclusions which yields that x* is 
a solution of the primal inclusion problem (|6j) and y* is a solution of the dual inclusion 
problem 0. 

Remark 6 As we have already mentioned, the system of inclusions ^ is solvable if and 
only if the primal inclusion problem ^ is solvable, which is further equivalent to solvability 
of the dual inclusion problem ([7]). Let us also notice that from the point of view of the 
numerical implementation Algorithm [4] has the drawback to ask for the calculation of the 
inverses of Idx +t 2 K*K and Idy +t 2 KK* . This task can be in general very hard, but 
it becomes very simple when K is, for instance, orthogonal, like it happens for the linear 
transformations to which orthogonal wavelets give rise and which play an important role 
in signal processing. 
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3 Zeros of sums of compositions of monotone operators with 
linear continuous operators 

In this section we provide via the primal-dual scheme Algorithm [T] an algorithm for solving 
the inclusion problem 

k 

find x £ X such that E ^UiK* BiKiX, (24) 

i=i 

where X and Yi are real Hilbert spaces, Bi : Yj zz£ Yj are maximally monotone operators, 
Ki : X — y Yi are linear and continuous operators for i = 1, k and u)i E (0, 1], i = 1, A;, 
are real numbers fulfilling X^i=i w « = !• The dual inclusion problem of (24) reads 

find y = (yi, ...,y fc ) e^x ... x Y fc such that ^^iVi = and (^(BiKi)' 1 fa) ^ 01 

i=l i=l 

(25) 

Following the product space approach from (8| (see also (5l) we show that this primal-dual 
pair can be reduced to a primal-dual pair of inclusion problems of the form (|6])-([7])- 

Consider the real Hilbert space H := X k endowed with the inner product (x,u)h = 
Eiti w t(a;i,itt)j for x = (xi)i<i<k,u = (ui)i<i< k € H, where {-,-)x denotes the inner 
product on X. Further, let Y := Y\ x ... X Y k be the real Hilbert space endowed with 
the inner product (y,z) Y ■= Yn^UiiVi, ^i)Y % for y = (yi)i<i<k,z = (zi)i<i<k £ Y, where 
{-,')Yi denotes the inner product on Yj,i = l,...,k. We define A : H z4 H, A := Ny, 
where V = {(x,...,x) £ H : x £ X}, B : Y =3 Y, B(y 1 ,...,y k ) = (B x y x , ...,B k y k ), and 
K : H — > Y, K(x\, Xk) = (K%xi, ...,K k x k ). Obviously, the adjoint operator of K 
is K* : Y -> H, K*( yi ,...,y k ) = {K* x y u K*y k ), for (yi,...,y fc ) E Y. Further, let be 
j : X H, j(x) = (x,...,x). 

The operators A and B are maximally monotone and 



x solves (|4j) if and only if (0, ...,0) E A(j{x)) + K*BK(j(x)), 

while 

y = (yi, ...,y fc ) solves @ if and only if (0, .., 0) E S^y - ^-^-iT^y. 
Applying Algorithm [T] to the inclusion problem 

find (xi, x k ) E H such that E ^4(xi, Xfc) + K* BK(x±, Xfc) (26) 
with starting point (x°, x°, y°, y°) E Xx...xX x Yj x ... x Y&, constants <r, r > 

s V ' 

fc 

and (a;^, := {x , ...,x°) yields for any n > the following iterations: 



i<j<fc 



(xr +1 )i<,< fc 



)Ki<k 



(y? +1 )i<i<k ■= J* B -i {(Vi)i<i<k + °K(x? 

" n+1 " ~JtA ((*?)i<i<* - rK* (y2 +1 )x<i< k ) 
2(x r l +1 ) 1 < i <k - (a?)i< f < fc . 

According to {.8), for the occurring resolvents we have that J t a{u\, -.-,u k ) = j(Yli=i WUi) 
for (u!,...,u k ) E H and J aB -i(zt, z k ) = {J B -\Z\, J B ~iz k ) for (zi,...,z fc ) E Y. This 

1 A; 
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x n k and x™ +1 



which shows 



means that for any n > 1 it holds x™ = .. 
that there is no loss in the generality of the algorithm when assuming that the first k 
components of the starting point coincide. Notice that a solution (xi, ...,Xk) of (26) must 
belong to domA, thus x\ = ... = x^. We obtain the following algorithm: 

Algorithm 7 

Initialization: Choose a, r > such that 0"r^i=i ll-^-ill 2 < 1 anc ^ 



For n > set: 



(x°,J/?,...,yg) G X x Yi x ... x Y fe . Set x 



o ._ 



y? +1 



2x 



n+1 



The convergence of Algorithm [7] is stated by the following result which is a consequence 
of Theorem [2j 



Theorem 8 Assume that X and Yi,i = 1, k, are finite- dimensional spaces and (24) is 



solvable. Then ( |25[ ) is ako solvable and the sequences (x ra ) n >o and (y™, ...,y]^) n >o generated 
in Algorithm^converge to a solution of the primal inclusion problem (24) and to a solution 
of the dual inclusion problem (25), respectively. 

Remark 9 Since \\K\\ 2 < Yli=i ll^«ll 2 ; the inequality o~T^2i=i ll-^«l| 2 < 1 i* 1 Algorithm [7] 
is considered in order to ensure that err 

ll^ll 2 < 1- 

When particularizing the above framework to the case when Yi = X and K{ = Idx for 
i = 1, .., k, the primal-dual pair of inclusion problems (24)-(25) become 



find i£l such that G 



i=i 



and 



find y = (yi,...,yk) £ X x ... x X such that ^^Wjyj 

i=l 



Oand p|Sri( y ^ 



(27) 



(28) 



i=l 



respectively. In this situation H = Y, K = Id#, ||-ftT|| = 1 and 



x solves (27) if and only if (0, ...,0) G A(j(x)) + B(j(x)), 



while 



{yi,...,y k ) solves pb if and only if (0, 0) G B~ l (y) - A~ l (-y). 



Algorithm [7] yields in this particular case the following iterative scheme: 



Algorithm 10 
Initialization: 



Choose a, r > such that ar < 1 and 

(x°,y°,...,y°) G X x ... x X. Set x° := . 



For n > set: y. 



n+1 



r n+l 
=-n+l 



J aB7 i(y? + k oS n ),i = 
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The convergence of Algorithm 10 follows via Theorem [8) 



Theorem 11 Assume that X is a finite- dimensional space and (27) is solvable. Then 
(28) is also solvable and the sequences (x") n >o and (y™, ...,y?) n >o generated in Algorithm 



10 converge to a solution of the primal inclusion problem (27) and to a solution of the 



dual inclusion problem (28), respectively. 



In the last part of this section we provide a second algorithm which solves (27) and 



( 28 ) which starts from the premise that by changing the roles of A and B one has 



x solves (27) if and only if (0, ...,0) G B(j(x)) + A(J(x)), 



while 



y=(yi,...,Vk) solves © if and only if (0, 0) G A~ l (—y) — B~ l (y). 



Applying Algorithm [T] to the inclusion problem 

find (xi, ...,Xk) 6 H such that £ B(x±, ...,Xk) + A(x\, x^) 



with starting point (x\ "'' 



Li 



, y^ 1 , y°) G X x ... x X, constants <r, r > and(x5, x^) 



2A; 



(x^, yields for any n > the following iterations: 



(yf + )\<%<k = JaA-i {{y?)i<i<k + o-(^r)i<i<fc 

W +1 )i<i<k = JrB((x2)i<i<k ~ r(yr +1 )i< i < fc ) 
(x? +1 )i<i<k = 2(x" +1 ) 1 < i < fc - (x?)i<i< fc . 

Noticing that J CTj 4-i = J^^" 1 = —°~Jcr-iN v a Wjj (cf. 5 Proposition 23.18]) and 

J (T -ij Vv ('Ui, ...,«fc) = JN v (ui,...,u k ) = j(X^LiWiiii) for (it 1} ...,itjfc) G il (cf. jsj relation 
(3.27)]) and by making for any n > the change of variables yf := — y™ for i = 1, A;, 
we obtain the following iterative scheme: 

Algorithm 12 

Initialization; Choose a, r > such that o~t < 1 and 

(x°,...,x°,y 1 °,...,y°)G Xx...xX . Set (x°,...,x°) := (*?,.., *°). 

2fc 



For n > set: y 



n+1 

i 

„n+l 



x? +1 



yf - < + EJ=1 "jltf + E -=i wjs?, i = 1, 
J TBi (xf + Ty? +1 ),i = l,...,k 



k 



2x 



n+1 



1,...,* 



Theorem 13 Assume that X is finite dimensional and (27) is solvable. Then (|28j) is also 
solvable and for all i = 1, A i/ie sequence (x") n >o generated in Algorithm\lS\ converges 
to a solution of (27) and i/ie sequence (y™ , y£) n >o generated by the same algorithm 
converges to a solution of (28). 
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4 Solving convex optimization problems via the primal-dual 
algorithm 



The aim of this section is to employ the iterative methods investigated above for solving 
several classes of unconstrained convex optimization problems. To this end we consider 
first for the real Hilbert spaces X and Y the proper, convex and lower semicontinuous 
functions / : X — )■ M and g : Y — > M and the linear and continuous operator K : X — >• Y 
the optimization problem 

inf {f(x) + g(Kx)} (29) 

ISA 

along with its Fenchel dual problem (see (o][d] [14|[24] ) 

SU p{-r(-K*y)-g*(y)}. (30) 

yeY 

For this primal-dual pair weak duality always holds, i.e., the optimal objective value of the 
primal problem is greater than or equal to the optimal objective value of the dual problem. 
In order to guarantee strong duality, i.e., the situation when the optimal objective values 
of the two problems coincide and the dual problem has an optimal solution one needs to 
ask for the fulfillment of a so-called qualification condition. Some of the most popular 
inferiority- type qualification conditions are (see, for instance, [5]~[7| [l4|[2T[[24] ) : 

(QCi) | 3x' G dom/ n K~ l (dom g) such that g is continuous at Kx', 
{QC 2 ) | G int(dom g - K(dom /)) 

and 

(QC3) I G sqri(dom# - K(domf)). 

We notice that (QC±) (QC2) => (QC3), these implications being in general strict, and 
refer the reader to the works cited above and the references therein for other qualification 
conditions in convex optimization. 

Algorithm [T] written for A := df and B := dg yields the following iterative scheme: 

Algorithm 14 

Initialization: Choose a,T > such that crr||ET|| 2 < 1 and (x°,y°) G X x Y. 
Set x° := x°. 
For n > set: y n+1 := prox CTff «(y n + aKx 11 ) 

x n+i ._ p r0 x T ^(x n - rK*y n+1 ) 

7pn+l ._ 2x n +l x n 



We have the following convergence result. 



Theorem 15 Assume that the primal problem (29) has an optimal solution x and one 
of the qualification conditions {QCi), i = 1,2,3, is fulfilled. Let (x n ,x n ,y n ) n >o be the 
sequence generated by Algorithm^I^ The following statements are true: 

(i) There exists y G Y, an optimal solution of the dual problem (30), the optimal 
objective values of the two optimization problems coincide and (x, y) is a solution of the 
system of inclusions 

Kx G dg*{y) and - K*y G df{x). (31) 



13 



(ii) For any n > it holds 
lUii ^||2 



2r 



+ (l-ar||^|| 2 )^- 



2ct 



< 



|z° — x\ l2 



2t 



+ 



\y° -y\\ 2 

2a 



(32) 



thus the sequence (x n ,y n ) n >o is bounded. 

(Hi) If X andY are finite- dimensional, then (x n ) n >o converges to an optimal solution 
of (29) and (y n ) n >o converges to an optimal solution of (30). 



Remark 16 (i) Statement (i) of the above theorem is well-known in the literature, (31) 
being nothing else than the system of optimality conditions for the primal-dual pair ( 29 )- 
(30) (see, for instance, (6, 14, 24|), while the other two statements follow from Theorem 

El 



(ii) The existence of optimal solutions of the primal problem ( 29 ) is guaranteed if, for 
instance, / is coercive and g is bounded below. Indeed, under these circumstances, the 
objective function of (29) is coercive and the statement follows via 24, Theorem 2.5.1(h)] 
(see, also, (Bl Proposition 15.7]). On the other hand, when / is strongly convex, then 
f+goK is strongly convex, too, thus (29) has an unique optimal solution (cf. [5j Corollary 
11.16]). 

(hi) We rediscovered above the iterative scheme and the convergence statement from [9] 
as a particular instance of the general results furnished in Section [2j 

For X and Y{ real Hilbert spaces, gi : Y{ — > M proper, convex and lower semicontinuous 
functions, K{ : X — > Y( linear and continuous operators and oji G (0,1], i = 1, real 
numbers fulfilling X]f=i w i = 1 consider the optimization problem 



inf } \oJigi{Kjx) 



(33) 



i=l 



and its Fenchel-type dual problem 



sup ^2~ujig*(yi)- 

ViEYi,i=l,...,k i=1 



(34) 



For the primal-dual pair (33)-(34) strong duality holds whenever one of the following 
qualification conditions is fulfilled (see, for instance, |6,8,24|): 

(QCf ) I 3x' £ Hi=i -Kj -1 (dom <?j) such that gi is continuous at Kix' , i = 1, k, 



and 



{QCf) 



G int (nLidomgj - {{Kix, K k x) : x G X}^j 



(QCf) | 0G sqri( Ui=l dom ff i-{ (K lX , K k x) 



Again, (QCf) =^ (QCf) => (QCf), the implications being in general strict. By taking 
Bi := dgi, i = 1, k, Algorithm [7] yields the following iterative scheme: 
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Algorithm 17 

Initialization: Choose a, r > such that ct^*L 1 ||-fQ|| 2 < 1 and 



(x°,y1, ...,yl) G X x Y x x ... x Y fc . Set x° := x° 



For n > set: y. 



n+l 
i 

5 n+l 



P rox CT g* {Vi + (TK i X n ),i = 1, 

22; n + l _ j; n 



-* n+l 



The convergence of Algorithm 17 is stated by the following result which is a conse- 
quence of Theorem [81 



Theorem 18 Assume that the primal problem (33) has an optimal solution x and one of 
the qualification conditions (QCf'), i = 1,2,3, is fulfilled. The following statements are 
true: 

(i) There exists (yi, ■ ■■,yk) £ Yi X ...xY&, an optimal solution of the dual problem ( |34[ ) ; 
the optimal objective values of the two optimization problems coincide and (x, y±, y^) is 
a solution of the system of inclusions 

k 

KiX G dg*(yi),i = 1, and ^ujiK*yi = 0. (35) 

i=l 

(ii) If X and Y are finite-dimensional, then the sequences (x n ) n >o and (yf, ■■■,y k l ) n >o 
generated in Algorithm 11 converge to an optimal solution of (33) and (34), respectively. 



problems (33) and (34) become 



Considering, finally, the particular case when Y\ = X and Ki = Idx, i = k, the 

fc 



inf J>&( a 



(36) 



and, respectively, 



sup ^-UigHVi 

yi£X,i=l,...,k i=1 



(37) 



The qualification conditions (QC| 



1, 2, 3, looks in this case like: 



(QC\ d ) | 3x' G Di=i dom gi such that ^ is continuous at x', i = 1, fc, 



(QCf) G hit ntl dom S* - {(*, x) : x G X} 



and, respectively, 



G sqri ( rii=i dom gi — {(x, x) : x G X} 



By particularizing Algorithm 17 we obtain: 

Algorithm 19 

Initialization: Choose a, r > such that o~t < 1 and 



.0 .,0 



jj?,..,j;°)g lx. . .xl . Set 



For n > set: y 



n+l 



fc+l 



r n+l 
prn+l 



P 10 *^*^ + crx n ),i = l,...,k 



15 



while Algorithm 12 gives rise to the following iterative scheme: 



Algorithm 20 

Initialization; Choose a, r > such that o~t < 1 and 

(x° 1 ,...,x k ,y 1 ,...,y° k )e X x . . . x X . Set (x?,...,xg) := (x° , 

2fc 



For n > set: y 



n+1 



x? +1 



prox Tft (x™ + ry™ +1 ),i = 1, k 
2x? +1 -x?,i = 

We have the following convergence theorem. 



jX],i 



1. 



1,...,* 



Theorem 21 Assume that the primal problem (36) /ias an optimal solution x and one of 
the qualification conditions (QC\ d ), 
true: 



1,2,3, is fulfilled. The following statements are 



(i) There exists (y±, yk) £ X x ... x X , an optimal solution of the dual problem ( |37[ ), 
the optimal objective values of the two optimization problems coincide and (x, yi, ...,yfc) is 
a solution of the system of inclusions 



x G dg*(yi 



1, k, and y^Ujyj 
i=l 



0. 



(38) 



(mJ //AT is finite- dimensional, then the sequences (x n ) n >o and (y", ■■■,y k ) n >o generated 
in Algorithm 19 converge to an optimal solution of (pBl and (37), respectively. 



(Hi) If X is finite- dimensional, then for all i = 1, k the sequence (x™) n >o generated 



in Algorithm 20 converges to an optimal solution of (36) and the sequence (y", ■■■yu) n >0 
generated by the same algorithm converges to an optimal solution of (37). 



Remark 22 One can notice that Theorem 12 1 1 remains valid even under a weaker condition 
than in (QC\ d ), namely by assuming that there exists x 1 G n* =1 doings such that k— 1 of 
the functions gi,i = 1, ...,k, are continuous at x' (see |6l Remark 2.5]). 



5 Numerical experiments 

In this section we present numerical experiments involving the primal-dual algorithm and 
some of its variants when solving some nondifferentiable convex optimization problems 
originating in image processing and in location theory. 



5.1 Image deblurring and denoising 

For a given matrix A £ l mxm describing a blur operator and a given vector b G M. m 
representing the blurred and noisy image the task that we considered was to estimate the 
unknown original image x* G M. m fulfilling 

Ax = b. 

With this respect we dealt with the regularized least squares problems 
(P 2 ) inf {aHxI^ + HAx-611 2 ) 
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and 



(P 3 



inf 



A 



+ \\Ax 



b\\ 2 + 5 s (x)}, 



where S C M m is an m-dimensional cube representing the range of the pixels and A > 
the regularization parameter. One of our aims was to show that in some concrete cases 
the quality of the recovered image via classical l\ regularized problem is by far not as good 
as the one recovered when regularizing with A|| • ||i + 8g. We solved problem (P2) by using 



Algorithm 14 and problem (P3) by using Algorithm 17 and showed the benefits of having 



the first one extended to problems having in their objective the sum of more than two 
functions. 

We concretely looked at the 272 x 329 blobs test image, which is part of the image 
processing toolbox in Matlab. We scaled the pixels to the interval [0, 1] and vectorized the 
image, obtaining a vector of dimension m = 272 x 329 = 89488. Further, by making use 
of the Matlab functions imf ilter and fspecial, we blurred the image as follows: 



H=fspecial ( ' gaussian ' ,9 ,4) ; % gaussian blur of size 9 times 9 

% and standard deviation 4 
B=imfilter (X,H, ' conv ',' symmetric ') ; % B=observed blurred image 

% X=original image 



In row 1 the function fspecial returns a rotationally symmetric Gaussian lowpass filter 
of size 9x9 with standard deviation 4. The entries of H are nonnegative and their sum 
adds up to 1. In row 3 the function imf ilter convolves the filter H with the image X and 
outputs the blurred image B. The boundary option "symmetric" corresponds to reflexive 
boundary conditions. 

Thanks to the rotationally symmetric filter H, the linear operator A £ ]^ mxm given 
by the Matlab function imf ilter is symmetric, too. By making use of the real spectral 

1. After adding a zero- mean white Gaussian 



\A\ 



decomposition of A it shows that 

noise with standard deviation 10~ 3 , we obtained the blurred and noisy image b 6 
which is shown in Figure [T] We solved the problem (P2) by applying Algorithm 14 for 



original 



blurred and noisy 




o 




Figure 1: The 272 x 329 blobs test image 



A = 2e - 6, a = 0.01, r = 9.99, / : 



Fill, 9 ■ 



k, g(y) = \\y-b\\- 
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and K = A. Since g*(y) = \\\y\\ 2 + (y, b) for y G M m , for all z G R m it holds 
prox^, (z) = argminjf ||y|| 2 + a(y,b) + ±\\z - y\\ 2 } = ^{z - ab), 



while 



prox r/ (z)j = argmin{TA|y|j + \\zi - yi\ 2 } = max{|zj| - rA, 0} sgn(zj) Vi = l,...,m. 



We solved the problem (P3) by applying Algorithm 17 for k = 3, oj{ = ^ , z = 1,2,3, 
A = 2e - 6, a = 0.05, r = 6.66, gi : R m -» K, gi (x) = X\\x\\i, K x = H R », 52 : M m -> M, 
92(y) = lb " &II 2 , ^2 = ^ and 53 : K m -> I,^) = 8g(x),K 3 = Id tt m. For all z G M m it 
holds 

P rox CT9l *( z )i = z i - cproxAn ^ (^)j = ^ - max{|zj| - A,0}sgn(zj) Mi = 1, ...,m, 

P rox <rg 2 *(^) = ^{z-crb) 

and 

Prox^^z) = z - crprox CT5g = z - aP s (±2) . 




PDlo = 2.4973 



PD? 00 = 0.4061 



PDf 50 = 0.1870 




Figure 2: Iterations 50, 100 and 150 for solving (P2) via Algorithm 14 and (P3) via 



Algorithm 17 



The top line of Figure [2] shows the iterations 50, 100 and 150 of Algorithm 14 for solving 
(P2); while the bottom line of it shows the iterations 50, 100 and 150 of Algorithm 17 for 
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solving (P3), for each of them the value of the objective function at the respective iterate 
being provided. All in all the quality of the recovered image by solving (P3) significantly 
outperformed the one of the image recovered by solving the classical l\ regularized least 
squares problem (P2). Moreover, in the images recovered by solving (P2) some artefacts 
could be identified. The gap between the quality of the recovered images is emphasized 
also by the improvement in signal-to-noise ratio (ISNR), which is defined as 



ISNR(n) = 10 log 



x ■ 



10 



\x — x'- 



where x, b and x n denote the original, observed and estimated image at iteration n, 
respectively. Figure [3] shows the evolution of the ISNR values when solving {P2) and (-P3). 
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A||.|| 1 + M(.)-6|| 2 

A||- h + \\A(-) - b\\ 2 + 5 mm (-) 



75 

Iterations 



150 



Figure 3: Improvement in signal-to- noise ratio (ISNR) 



5.2 The Fermat- Weber problem 

The second application of the primal-dual algorithm presented in this paper is with respect 
to the solving of the Fermat- Weber problem, which concerns the finding of a new facility 
in order to minimize the sum of weighted distances to a set of fixed points. We considered 
the nondifferentiable convex optimization problem 



{Pfw) inf ly^Xi 



where Cj G M. m are given points and Aj > are given weights for i = 1, k. We solved the 



optimization problem (Pfw) by using Algorithm 19 for = | and gi : W 71 — > M.,gi(x) 



Xi\\x — Ci\\,i = 1, k. With this respect we used that for i = 1, k it holds 

m = { ^ if i yll " Aj ' Vy G M. m 
^ 1 +00, otherwise, 
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and, from here, when a > 0, that 



prox, 



A 



z - aa, 

z—aci 



if \\z — aci\\ < Xi, 
otherwise 



Vz e 



We investigated the functionality of the algorithm on two prominent sets of points and 
weights, often considered in the literature when analyzing the performances of iterative 
schemes for the Fermat- Weber problem. 

In a first instance we considered for k = 4 the points in the plane and the weights 

a = (59, 0), c 2 = (20, 0), c 3 = (-20, 48), c 4 = (-20, -48) and Ai = A 2 = 5, A 3 = A 4 = 13, 

(39) 

respectively. The optimal location point is x = (0,0), however, the classical Weiszfeld 
algorithm (see 17,23 ) with starting point x° = (44,0) breaks down in (20,0). On the 
other hand, Algorithm 19 with a = 0.13, r = 1.4 and = (0,0), k = 1,...,4, achieved a 



point which is optimal up to three decimal points after 30 iterations. Figure [4] shows the 
progression of the iterations, while PD n provides the value of the objective function at 
iteration n. Recently an approach for solving the Fermat- Weber problem was proposed by 



PD = 2275.0 



PD 3 = 1891.8 




-20 -10 10 20 30 40 50 

PD 5 = 1777.94 



20 -10 10 20 30 40 50 60 

PD 30 = 1747.0 




-20 -10 10 20 30 40 50 60 



-20 -10 10 20 30 40 50 60 



Figure 4: The progression of iterations of Algorithm 19 when solving the Fermat- Weber 



problem for points and weights given by (39) 



Goldfarb and Ma in [16] , which assumes the approximation of each of the functions in the 
objective by a convex and differentiable function with Lipschitz-continuous gradient. The 
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optimization problem which this smoothing method yields is solved in [TBI by the classical 
gradient method (Grad) and by a variant of Nesterov's accelerated gradient method (Nest) 
(see uM) and by a fast multiple-splitting algorithm (FaMSA-s) introduced in this paper. 
We applied the smoothing approach in connection with these algorithms to the example 



considered in (39) with smoothness parameter p equal to 10~ 3 (chosen also in 16 ) and 
step sizes r = 0.1, r = 0.01 and r = 0.001. We stopped the three algorithms when 
achieving an iterate x n such that \\x n — x\\ < 10~ 3 and obtained in all cases the lowest 
number of iterations for r = 0.1. A point which is optimal up to three decimal points was 
obtained for Nest after 308 iterations, for Grad after 175 iterations and for FaMSA-s after 



54 iterations, none of these iterative schemes attaining the performance of Algorithm 19 
For the second example of the Fermat- Weber problem we considered in case k = 5 the 
points in the plane and the weights (see [12|) 



ci = (0, 0), c 2 = (1, 0), c 3 = (0, 1), c 4 = (1, 1), c 5 = (100, 100) 
and Ai = A2 = A3 = A4 = 1, A5 = 4, 

respectively. The optimal location point is x = (100, 100) and, by choosing the relative 
center of gravity x° = (50.25, 50.25) as starting point, we found out that not only the clas- 
sical Weiszfeld algorithm, but also the approach from [16] described above in connection 
to each of the methods Grad, Nest and FaMSA-s did not achieve a point which is optimal 



up to three decimal points after millions of iterations. On the other hand, Algorithm 19 
with a = 0.0001, r = 9999 and y® = (0, 0), k = 1, 5, achieved a point which is optimal 
up to three decimal points after 478 iterations. This example is more than illustrative 



for the performance of the primal-dual Algorithm 19 in comparison to some classical and 
recent algorithms designed for the Fermat- Weber problem. Figure [5] shows the progression 
of the iterations, while PD n provides the value of the objective function at iteration n. 



6 Conclusions 

In this paper we motivate and formulate a primal-dual algorithm which solves both the 
problem of finding the zeros of the sum of a maximally monotone operator with the 
composition of another maximally monotone operator with a linear continuous operator 
and its Attouch-Thera-type dual inclusion problem in Hilbert spaces. We also investigate 
the convergence of the provided iterative scheme and show how one can derive from it a 
splitting algorithm for finding the zeros of the sum of compositions of maximally monotone 
operators with linear continuous operators. As particular instances of the general schemes 
algorithms for solving several classes of nondifferentiable convex optimization problems 
are introduced. Among them one can rediscover the primal-dual algorithm from |9J for 
solving the problem which assumes the minimization of the sum of a proper, convex 
and lower semicontinuous function with the composition of another proper, convex and 
lower semicontinuous function with a linear continuous operator. The performances of the 
provided algorithm are emphasized in the context of some applications in image deblurring 
and denoising and in location theory. 

Acknowledgements. The authors are thankful to Christopher Hendrich for the imple- 
mentation of the numerical schemes to which comparisons of the primal-dual algorithm 
were made. 
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Figure 5: The progression of iterations of Algorithm 19 when solving the Fermat- Weber 



problem for points and weights given by (40) 
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